Statistics of mixing in three-dimensional Rayleigh— Taylor 
turbulence at low Atwood number and Prandtl number one 



G. Boffetta^, A. Mazzino^, S. Musacchio^ and L. Vozella^ 

^ Dipartimento di Fisica Generate and INFN, 
Universitd di Torino, via P.Giuria 1, 10125 Torino (Italy) 
and CNR-ISAC, Sezione di Torino, 
corso Fiume 4, 10133 Torino (Italy) 
Dipartimento di Fisica, Universitd di Genova, 
INFN and GNISM, via Dodecaneso 33, I6I46 Genova (Italy) and 
3 GNRS UMR 6621, Lab. J. A. Dieudonnne, 

Universite de Nice Sophia- Antipolis, 
Pare Valrose, 06108 Nice Gedex 02 (France) 
(Dated: June 23, 2010) 

Abstract 

Three-dimensional miscible Rayleigh-Taylor (RT) turbulence at small Atwood number and at 
Prandtl number one is investigated by means of high resolution direct numerical simulations of the 
Boussinesq equations. RT turbulence is a paradigmatic time-dependent turbulent system in which 
the integral scale grows in time following the evolution of the mixing region. In order to fully 
characterize the statistical properties of the flow, both temporal and spatial behavior of relevant 
statistical indicators have been analyzed. 

Scaling of both global quantities [e.g., Rayleigh, Nusselt and Reynolds numbers) and scale de- 
pendent observables built in terms of velocity and temperature fluctuations are considered. We 
extend the mean-field analysis for velocity and temperature fluctuations to take into account in- 
termittency, both in time and space domains. We show that the resulting scaling exponents are 
compatible with those of classical Navier-Stokes turbulence advecting a passive scalar at compara- 
ble Reynolds number. Our results support the scenario of universality of turbulence with respect 
to both the injection mechanism and the geometry of the flow. 



I. INTRODUCTION 



The Rayleigh-Taylor (RT) instability is a well-known fluid-mixing mechanism originating 
at the interface between a light fluid accelerated into an heavy fluid. It was flrst described 
by Rayleigh |l| for incompressible fluid under gravity and later generalized to all accelerated 
fluid by Taylor fl. 

RT instability plays a crucial role in many flelds of science and technology. In particular, 
in gravitational fusion it has been recognized as the dominant acceleration mechanism for 
thermonuclear reactions in type-la supernovae 0, 4|. The efficiency of inertial confinement 
fusion depends dramatically on the ability to suppress RT instability on the interface between 
the fuel and the pusher shell {5,6]. 

In a late stage, RT instability develops into the so-called RT turbulence in which a layer 
of mixed fiuid grows in time increasing the kinetic energy of the fiow at the expenses of 
the potential energy. This process finds applications in many fields, e.g. atmospheric and 
oceanic buoyancy driven mixing. Despite the great importance and long history of RT 
turbulence, a consistent phenomenological theory has been proposed only recently In 
three dimensions, this theory predicts a Kolmogorov-like scenario, with a quasi- stationary 
energy cascade in the mixing layer. The prediction is based on the Kolmogorov-Obukhov 
picture of turbulence in which density fiuctuations are transported passively in the cascade 
and kinetic-energy fiux is scale independent js]. Quasi-stationarity is a consequence of Kol- 
mogorov scaling of characteristic times associated to turbulent eddies: large-scales grow 
driven from potential energy, while small-scale structures, fed by the turbulent cascade, 
follow adiabatically large-scale growth. These theoretical predictions have been partially 
confirmed by recent numerical studies js], ojjlll. Other alternative phenomenological ap- 
proaches (see e.g. 12|) does not necessarily lead to the Kolmogorov scaling for the energy 
spectra. 

In this Paper we carry out an analysis of the scaling behavior of relevant observables with 
the aim of deepening our previous investigation Indeed, our aim is to make a careful 
investigation of the time evolution of global observables and of spatial/temporal scaling and 
intermittency. Moreover we push the analogy of RT turbulence with usual Navier-Stokes 
(NS) turbulence much further. We show that small-scale velocity and temperature fiuctua- 
tions develop intermittent distributions with structure-function scaling exponents consistent 
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with NS turbulence advecting a passive scalar. 

This Paper is organized as follows. In Sec. [Ill we formulate the problem and outline the 
phenomenology. After providing a description of the numerical setup in Sec. Illlt we describe 
our results in the subsequent Sections. Sec. |IV]is devoted to the investigation of the temporal 
evolution of global quantities. In Sec. |V] we focus on the statistics at small scales. Finally, 
the Conclusions are provided by summarizing the main results. 

II. EQUATION OF MOTION AND PHENOMENOLOGY 

We consider the three-dimensional Boussinesq equations for an incompressible velocity 
field (V ■ V = 0), 

(9tv + V • Vv = -Vp + z/Av - (3gT (1) 
dtT + v-VT = kAT (2) 

T(x, t) being the temperature field, proportional to the density via the thermal expansion 
coefficient (3 a.s p = po[l — f3(T — To)] (po and Tq are reference values), u is the kinematic 
viscosity, k, the molecular diffusivity and g = (0, 0, —g) the gravitational acceleration. 

At time t = the system is at rest with cooler (heavier, density P2) fluid placed above the 
hotter (lighter, density pi) one. This corresponds to v(x, 0) = (0, 0, 0) and to a step function 
for the initial temperature profile: T(x, 0) = — (^^o/2)sgn(2;) where 6q is the temperature 
jump which fixes the Atwood number A = (p2 — Pi)/(P2 + P2) = (l/2)/3^o- The development 
of the instability leads to a mixing zone of width h which starts from the plane z = 
and is dimensionally expected to grow in time according to h{t) = aAgt^ (where a is a 
dimensionless constant to be determined) which implies the relation Vrms — Agt for typical 
velocity fluctuations (root mean square velocity) inside the mixing zone. 

The convective state is characterized by the turbulent heat flux and energy transfer as a 
function of mean temperature gradient. In terms of dimensionless variables these quantities 
are represented respectively by the Nusselt number Nu = 1 + {wT)h/ {k^Oq) {w being the 
vertical velocity) and the Reynolds number Re = Vrmsh/v as a function of the Rayleigh 
number Ra = figO^h^ / {vk) and the Prandtl number Pr = u/n. Here and in the following 
(...) denotes spatial average inside the turbulent mixing zone, while the overbar indicates 
the average over horizontal planes at flxed z. 
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One of the most important problems in thermal convection is to find the functional 
relation between the convective state characterized by Nu and Re and the parameter space 
defined by Ra and Pr 13|. The existence of an asymptotic regime at high i?a, with a 
simple power law dependence Nu ~ Ra^ and Re ~ Ra' , is still controversial in the case 
of Rayleigh-Benard convection, despite the number of experiments at very large Ra. Most 
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of a more complex behavior 



of the experiments have reported an exponent ^ ~ 0.3 

[l6 . IitI partially described by a phenomenological theory [l8|. However, many years ago, 
Kraichnan lol predicted an asymptotic exponent ^ = 1/2 (with logarithmic corrections) 
associated to the now called "ultimate state of thermal convection", while exponents ^ > 1/2 
are excluded by a rigorous upper bound Nu < (l/6)-Ra^/^ — 1 [20]. The ultimate state regime 
is expected to hold when thermal and kinetic boundary layers become irrelevant, and indeed 
has been observed in numerical simulations of thermal convection at moderate Ra when 
boundaries are removed 211, while no indication of ultimate state regime has been observed 

n 

in Rayleigh-Benard experiments |14| . 

The ultimate state exponent is formally derived from kinetic energy and temperature 
balance equations [l8|. In the present context of RT turbulence they can more easily be 
obtained from the temporal scaling of h and Vrms- Assuming that (wT) ~ Vrms^o, using the 
above definitions one estimates: 



Ra ~ {Ag)Hy{i^K), Re ~ (A^)V/z/ and Nu ~ (A^)V/k 



from which 



Nu ~ Pr'/'Ra'/' 



and Re ~ Pr-^/^Ra^'"^ 



(3) 



(4) 



For what concerns the small-scale statistics inside the mixing zone, the phenomenological 
theory j?! predicts for the 3D case an adiabatic Kolmogorov-Obukhov scenario with a time- 
dependent kinetic-energy fiux e ~ "'^ ih ^ ( R^aA'i 
and temperature fiuctuations are therefore expected to follow 



~ "^rms/^ — (/^fl'^o) ^- Spatial-temporal scaling of velocity 



5rT{t) ~ e-i/^ei/V^/^ ^ el''\Pg)-^'H-^/''r^''' 



(5) 
(6) 



where 5rv{t) = v{x + r,t) — v{x, t) is the velocity increment on a separation r (similarly for 
temperature) and ^ do^t~^ is the temperature- variance fiux. We remark that the above 
scaling is consistent with the assumption of the theory that temperature fiuctuations are 



Label 


N.,=Ny 


Nz I' = K, 




A 


256 


1024 9.5 X 10"^ 


103 


B 


512 


2048 4.8 X 10-^ 


196 


C 


1024 


1024 3.2 X 10"^ 


122 



TABLE I: Parameters of the simulations. N^, Ny, spatial resolution, viscosity, k thermal 



difFusivity, 



,y^l5/(i^e) Reynolds number evaluated at the end of the simulation. All 



dimensional quantities are made dimensionless using the vertical box size L^, the characteristic 
time T = (Lz/Ag)^/"^ and the temperature jump as reference units. 



passively transported at small scales (indeed using (jSIH]) the buoyancy term (3gT becomes 
subleadinff in (IT]) at small scales). This is the main difference with respect to the 2D case 

n p n 

in which temperature fluctuations force the turbulent flow at all scales |7|, |22| . 



III. NUMERICAL SETTING 



The Boussinesq equations ([HE]) are integrated by a standard 2/3-dealiased pseudospectral 
method on a three-dimensional periodic domain of square basis = Ly and aspect ratio 
L^/Lz = R with uniform grid spacing at different resolutions as shown in Table [T] In the 
following, all physical quantities are made dimensionless using the vertical scale Lz, the 
temperature jump 9q and the characteristic time r = (Lz/ AgY^'^ as fundamental units. 

Time evolution is obtained by a second-order Runge-Kutta scheme with explicit linear 
part. In all the runs, I3g = 2.0 and Pr = u/ k, = 1. Viscosity is sufficiently large to resolve 
small scales (kmaxV — 1-2 at final time, being r] = u^/^e~^^^ the Kolmogorov scale and 

RT instability is seeded by perturbing the initial condition with respect to the unstable 
step profile. Two different perturbations were implemented in order to check the indepen- 
dence of the turbulent state from initial conditions. In the first case the interface T = 
at 2 = is perturbed by a superposition of two-dimensional waves of small amplitude 
ho = O.OOALz in an isotropic range of wavenumbers 32 < /c < 64 (with k"^ = k"^ + ky) and 
random phases |23j. For the second set of simulations, we perturbed the initial condition 
by adding 10% of white noise to the value of T(x, 0) in a layer of width around z = 0. 
Figure [I] shows a snapshot of the temperature field in a cubic slice around z = m the 
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turbulent regime at time t = 2t for simulation B (see Table [T]). 




FIG. 1: Snapshot of temperature field for Ray leigh- Taylor simulation at t = 2r. White (black) 
regions corresponds to hot (cold) fluid. Parameters in Table [U run B. 

IV. EVOLUTION OF GLOBAL QUANTITIES 

Figure |2] displays the evolution of the total kinetic energy E = J (l/2)t;(x)2 rfx and total 
kinetic-energy dissipation e/, as a function of time. After the linear instability regime, at 
t ~ r the turbulent regime sets in with algebraic time dependence. Temporal evolution of 
the two quantities are easily obtained recalling that, being global quantities, an additional 
geometrical factor h{t) ~ t"^ due to the integration over the vertical direction has to be 
included. Therefore the predictions are E{t) ~ vl^^h ~ and ~ e/i ~ t^, as indeed 
observed at late times. We also plot in Fig. [2] the total potential-energy loss, defined as 
P(0) — P{t) with P{t) = —Pg J zT{x.) (ix which has the same temporal scaling of E{t) as 
it is obvious from energy balance: d{E + P)/dt = — e^. Notice that for this non-stationary 
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FIG. 2: Temporal growth of kinetic energy E (red circles), kinetic-energy dissipation (blue 
triangles) and potential-energy loss AP (pink squares) for run B. For clarity of the plot cl has 
been shifted by a factor 10. The two short lines represent the dimensional scaling E{t) ~ AP(t) ~ 
and eiit) ~ t^. Inset: ratio of the energy growth rate dE/dt and the flux ei- Data from run B. 

turbulence the energy balance does not fix the ratio between the energy growth rate dE/dt 
and the energy dissipation (and flux) e/,. In the turbulent regime, our simulations show 
an "equipartition" between large-scale energy growth and small-scale energy dissipation: 
dE/dt c:^ ei —{1/2) dP/dt. This amounts to saying that half of the power injected into 
the flow contributes to the growth of the large-scale flow, and half feeds the turbulent cascade 
(see inset of Fig. |2]). This result was found to be independent on the value of viscosity (the 



only adjustable parameter in the system) and is consistent with previous flndings 24 1. 

An interesting remark is that RT turbulence represents an instance of the general case 
of a turbulent flow adiabatically evolving under a time- dependent energy input density I{t) 
which forces the flow at the integral scale L(t) (concerning the problem of turbulent flow 



26( 1 and references therein). 



characterized by a time dependent forcing see, for example, 
Energy balance requires dS/dt = I{t) — e{t), where S is the kinetic energy density. Assuming 
a Kolmogorov spectrum for velocity fluctuations at scales smaller than the integral scale, 
one estimates £{t) ^ e^l^E^I"^ . Therefore, in situations characterized by an algebraic growth 
of the energy input density X(t) ~ t'^ a self-similar evolution of the energy spectrum can 
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be obtained only if e(t) ~ t"' and L{t) ~ t(3+7)/2 This is indeed realized in RT turbulence, 
where 7 = 1 and e ~ t, L(t) ~ t^. 

In the inset of Fig. |3] the growth of vertical and horizontal rms velocity {Wrms and Urms 
respectively), computed within the mixing layer, is shown. Both Urms and Wrms grow linearly 
in time, as expected, with the vertical velocity about twice the horizontal one, reflecting the 
anisotropy of the forcing due to gravity. It is interesting to observe that anisotropy decays 
at small scales, where almost complete isotropy is recovered, as shown in Fig. [3J The 
ratio of vertical to horizontal rms velocity reaches a value Wrms/urms — 1-8 at later times 
(corresponding to Rx ~ 200) while for the gradients we have {dzw)rms/{dxu)rms — l-O- 




20 40 60 80 100 120 140 160 180 200 



FIG. 3: Ratio of the vertical rms velocity Wrms to the horizontal rms velocity Urms (red open circles) 
and ratio of the vertical velocity gradient {dzw)rms to the horizontal velocity gradient {dxu)rms 
(blue filled circles) versus Reynolds number R\ indicating the recovery of isotropy at small scales. 
Inset: temporal evolution of horizontal rms velocity Urms (red open circles) and vertical rms velocity 
Wrms (blue filled circles). The black line represents linear scaling. Data from run B. 



The evolution of the mean temperature profile T(z, t) = l/(Lr.Ly) J T(x, t)dxdy is shown 



in Fig. m As observed in previous simulations [10|, lUl, l22|, |27| the mean profile is ap- 
proximately linear within the mixing layer (where therefore the system recovers statistical 
homogeneity). Nevertheless, statistical fluctuations of temperature in the mixing layer are 
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FIG. 4: Mean temperature profiles T{z, t) for a single realization of simulation B with diffused 
initial perturbation at times t = 1.4r, t = 2. Or, t = 2.6t and t = 3.2t. Lower and upper insets: 
profiles of the heat flux wT{z,t) and square vertical velocity w'^{z,t) at times t = 1.4r, t = 2. Or 
and t = 2.6r. 

relatively strong: at later time we find a fiat profile of fluctuations. Moreover their distri- 
bution is close to a Gaussian with a standard deviation (Tt{z) ~ 0.25^o (not shown here). 

In Fig. m we also plot the proflle of the heat flux wT{z, t) and the square vertical velocity 
w'^{z,t). Both vanish outside the mixing layer and inside show a similar shape not far from 
a parabola. Of course, the time behaviors of the heat-flux and of the square vertical velocity 
amplitude are different. Indeed, the former is expected to grow as oc t and the latter as oc t^. 

The mean temperature proflle deflnes the width of the mixing layer. Different deflnitions 
of the mixing width, /i, have been proposed on the basis of integral quantities or threshold 



values (see [28| for a discussion of the different methods). In the following we will use the 
simple deflnition based on a threshold value: T{±h/2) = s9o/2 where s < 1 represents the 
threshold. 

The evolution of the mixing width for s = 0.8 is shown in Fig. [51 After an initial 
stage {t < 0.3 r) in which the perturbation relaxes towards the most unstable direction, we 
observe a short exponential growth corresponding to the linear RT instability. At later times 
[t > 0.6 r) the similarity regime sets in and the dimensional t"^ law is observed. The nai've 
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FIG. 5: Evolution of the mixing-layer width h as a function of time t for simulation B computed 
from the profiles of Fig. H] with a threshold s = 0.8. The inset shows the compensation with 
dimensional prediction h/{Agt^) converging to a value ~ 0.036. 

compensation with Agt^ gives an asymptotic constant value h/{Agt^) ~ 0.036 for t > 3r 
and Re ~ 10^ (at which the mixing width is still below half box). For the calculation of a, 



more sophisticated analysis have been proposed recently J3|, |29|, |30|. using slightly different 
approaches (briefly, in j29| a similarity assumption and in jSOj a mass flux and energy balance 
argument). In both cases, the authors derive for the evolution of h{t) the equation 



= AaAgh 



(7) 



which has solution h{t) = aAgt^ + 2{aAhQYl'^t + Hq where Hq is the initial width introduced 
by the perturbation, a = h'^/{AAgh). The idea is to get rid of the subleading terms and 
extract the contribution at early time by using directly (j7]) and evaluating a = / (AAgh). 

The growth of the mixing layer width h{t), a geometrical quantity, is accompanied by the 
growth of the integral scale L{t), a dynamical quantity representing the typical size of the 
large-scale turbulent eddies. Following Ref. joj] we define L as the half width of the velocity 
correlation function f{L) = {vi{r)vi{r + L)) / (f^) = 1/2. In the turbulent regime the integral 
scale and the mixing length are linearly related (see Fig|6]). A linear fit gives L/h 1/17 and 
L/h ~ 1/42 for the integral scale based on the vertical and horizontal velocity component 
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respectively, in agreement with the results shown in jo] (of course, the precise values of the 
coefficients depend on the definition of h) . The anisotropy of the large scale fiow is refiected 
in the velocity correlation length: the integral scale based on horizontal velocity is smaller 
than the one based on vertical velocity. 




FIG. 6: Growth of the integral scale L based on the vertical velocity (red circles), and horizontal 
velocities (blue squares) as a function of the mixing layer h. Data from simulation B. 

We end this Section by discussing the behavior of the turbulent heat fiux, the energy 
transfer and the mean temperature gradient in terms of dimensionless variables (as dis- 
cussed in Sec [TTl): Nusselt, Reynolds and Rayleigh numbers, respectively. The temporal 
evolution of these numbers, shown in Fig. [TJ follows the dimensional predictions (|3]) for the 
temporal evolution of a (see Inset of Fig. |5]). The presence of the "ultimate state of thermal 
convection", in the restricted case Pr = 1, is also confirmed by our numerical results. Data 
obtained from simulations at various resolution (see Fig. [S]) are in close agreement with the 
"ultimate state" scalings @. 



V. SMALL-SCALE STATISTICS 



As already discussed in the introduction, the phenomenological theory predicts that, at 
small-scales, RT turbulence realizes an adiabatically evolving Kolmogorov-Obukhov scenario 
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FIG. 7: Temporal scaling of Nusselt number Nu = 1 + {wT)h/ [kOq) (blue triangles), Reynolds 
number Re = Vrmsh/v (black squares) and Rayleigh number Ra = iSgO^h^ /{uk) (red circles) for 
simulation B at Pr = 1. The lines are the temporal scaling predictions for Nu and Re and 
for i2a. 

of NS turbulence. Here adiabatic means that, because of the scaling laws, small scales have 
sufficient time to adapt to the variations of large scales, leading to a scale-independent 
energy flux. We remark that this is not the only possibility, as in two dimensions the 
phenomenology is substantially different. Unlike the 3D configuration, the 2D scenario is an 
example of active scalar problem. Indeed, the buoyancy effect is leading at both large and 
smaller scales. An adiabatic generalization of Bolgiano-Obukhov scalin g ha s been predicted 
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by means of mean field theory [3] and has been confirmed numerically 

Figure [9] shows the global energy fiux in spectral space at different times in the turbulent 
stage of the simulation. As discussed above, the flux grows in time following the increase 
of the input X{t) at large scales and at smaller ones, faster scales have time to adjust their 
intensities to generate a scale independent flux. 

If the analogy with NS turbulence is taken seriously, one can extend the dimensional 
predictions ([MS]) to include intermittency effects. Structure functions for velocity and tem- 



12 




FIG. 8: Nusselt (blue empty symbols) and Reynolds (black filled symbols) numbers as a function 
of Rayleigh number at Pr = 1. The "ultimate state" prediction (black and blue lines) expressed 
by ([H) are compared with numerical data obtained from simulations A (squares), B (circles) and 
C (triangles). 



perature fluctuations are therefore expected to follow 



S^{r,t) ^ ((M(t))n ^ ' (9) 

In ([8]) we introduce the longitudinal velocity differences SrVi\{t) = (v(x + r, t) — v(x, t)) -r/r 
and the increment r is made dimensionless with a characteristic large scale which, in the 
present setup, is proportional to the width of the mixing layer h{t), the only scale present in 
;he system. The two sets of scaling exponents (p and are known from both experiments 
3ll . l32| and numerical simulations 33| with good accuracy for moderate p. Mean-field 
prediction is Cp = Cj = p/3 while intermittency leads to a deviation with respect to this 
linear behavior. Kolmogorov's "4/5" law for third-order velocity implies the exact result = 
1, while temperature exponents are not fixed, apart for standard inequality requirements 
jSj. Both experiments and simulations give stronger intermittency in temperature than in 
velocity fluctuations, i. e. < Cp for large p. 
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FIG. 9: Spectral global kinetic energy flux n(/c) at times t = 2. At, t = 2.6t, t = 2.8t, t = 3. Or 
(from bottom to top) and temperature variance flux at t = 3. Or (inset) for simulation B. Kinetic 



energy flux is defined as n(A;) = — Re 'Uj(— k')(v • Vvi)(k') 



dk' where is the Fourier transform 



BJ. A similar definition holds for the temperature variance flux. 



We have computed velocity and temperature structure functions and spectra in our sim- 
ulations of RT turbulence. To overcome the inhomogeneity of the setup, velocity and tem- 
perature differences (at fixed time) are taken between points both belonging to the mixing 
layer as defined above. Isotropy is recovered by averaging the separation r over all direc- 
tions. Spectra are computed by Fourier-transforming velocity and temperature fields on 
two-dimensional horizontal planes and then averaging vertically over the mixing layer. 



A. Lower-order statistics 



In Fig. [TOT a). we plot kinetic-energy spectra at different times in the turbulent stage, 
compensated with the time dependent energy dissipation e^/^(t). In the intermediate range 
of wavenumbers, corresponding to inertial scales, the collapse is almost perfect. The evo- 
lution of the compensated spectra shows that the growth of the integral scale at small 
wavenumbers is in agreement with Fig. [61 Likewise temperature-variance spectra are con- 
sidered in Fig. [TOT b). Here, the spectra are compensated with both the time dependent 
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FIG. 10: (a) Kinetic-energy spectra compensated with e^/^ at times t = It (red crosses), t = 1.4r 
(green times), t = 1.8r (blue stars) and t = 3.8r (pink squares). Inset: kinetic-energy dissipation vs. 
time. The Hue represents the hnear growing of energy dissipation (see Sec. |Tlj). (b) Temperature- 
variance spectra compensated with e^^e^/^ at same times. Inset: temperature-variance dissipation 
vs. time. The hne is the dimensional prediction ~ (see Sec. [Tij). Data from simulation B. 
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FIG. 11: Third-order isotropic longitudinal velocity structure function 5*3 (r) computed at a late 
stage in the simulation (red circles) and mixed longitudinal velocity-temperature structure function 
'S'i,2(^) (blue triangles). The black line represents the linear scaling. Data from simulation B. 

temperature variance dissipation e^^(t) and the energy dissipation e^^^{t). The evolution of 

15 



the intermediate range of wavenumbers follows the dimensional prediction (|6]). 

Figure [TT] displays the third-order velocity structure function 6*3 (r), related to the energy 
flux by Kolmogorov's "4/5" law S'3(r) = — (4/5)er js]. We also plot the mixed velocity- 
temperature structure function 5'i_2(?") = {5rV\\ (5,.T)^) which is proportional to the (constant) 
flux of temperature fluctuations e-r according to Yaglom's law 5'i_2('") = — (4/3)eTT' |34|. 
Both the computed structure functions display a range of linear scaling, i.e. a constant 
flux, in the inertial range of scales 5 x 10"'^ ^ r / < 5 x 10~^. It is interesting to observe 
that the mixed structure function 5*1,2 1'") seems to have a range of scaling which extends to 
larger scales. This is probably due to the fact that at large-scale temperature fluctuations 
are dominated by unmixed plumes which have strong correlations with vertical velocity. 



Q. 



Q. 




FIG. 12: Structure-function scaling exponents for velocity increments ("p (circles) and temperature 
increments ("J (triangles) with absolute values. Red open symbols are obtained using ESS pro- 



cedure 35|] on the present simulation at time t = 3t, fixing the value of Cs = 1 and = 2/3. 
Errors represent fluctuations observed in difl^erent realizations of simulation B. Blue filled sym- 
bols are taken from a stationary NS simulation at Rx = 427 [33]. Black line is Kolmogorov 
non-intermittent scaling p/3. Insets: probability density function for velocity diff'erences 6rv{t) 
(upper) and temperature differences 5rT{t) at time t = 3r and scales r = O.OOSL^ (red circles) and 



r = O.O6L2 (green squares). Black lines represent a standard Gaussian. 
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B. Spatial/temporal intermittency 



Despite the clear scaling observable in Fig. [TTl it is very difficult to compute scaling ex- 
ponents directly from higher-order structure functions because of limited Reynolds number 
and statistics. Therefore, assuming a scaling region as in Fig. [TTl 'we can compute relative 
scaling exponents using the so-called Extended Self Similarity procedure 35|. This corre- 
sponds to consider the scaling of one structure function with respect to a reference one {e.g. 
Ss{r) for velocity statistics), and thus to measure a relative exponent {i.e. Cp/Cs)- 

Scaling exponents obtained in this way are shown in Fig. [121 Reference exponents for 
the ESS procedure are Cs = 1 and (2 = 2/3 (which is not an exact result). We see that 
both velocity and temperature scaling exponents deviate from the dimensional prediction of 
([MS]) {i.e. Cp = Cj = P/3) indicating intermittency in the inertial range. We also observe 
a stronger deviation for temperature exponents, which is consistent with what is known for 
the statistics of a passive scalar advected by a turbulent flow [s, [s^. 

The question regarding the universality of the set of scaling exponents with respect to the 
geometry and the large-scale forcing naturally arises. Several experimental and numerical 
investigations in three-dimensional turbulence support the universality scenario in which 
the set of velocity and passive-scalar scaling exponents are independent of the details of 
large-scale energy injection and geometry of the flow. Therefore, because we have seen 
that in 3D RT turbulence at small scales temperature becomes passively transported and 
isotropy is recovered, one is tempted to compare scaling exponents with those obtained in 
NS turbulence. As shown in Fig. [T21 the two sets of exponents coincide, within the error 
bars, with the exponents obtained from a standard NS simulation with passive scalar at 
comparable -Ra [33|. 

We remark that scaling exponents for passive scalar in NS turbulence are very sensitive to 
the fitting procedure. Strong temporal fluctuations have been observed in single realization 
37l | and dependence on the fitting region has been reported [33| . Indeed, different realizations 
of RT turbulence (starting with slightly different initial perturbations) lead to fluctuations 
of scaling exponents which account for the errorbars shown in Fig. [T^ 

Figure [121 also shows probability density functions for velocity and temperature fluctua- 
tions for two different scales. Both distributions are close to a Gaussian at large scale and 
develop wide tails at small scales, indicating the absence of self-similarity thus con&ming 
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the intermittency scenario. 
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FIG. 13: Time dependence of p-order velocity structure function Sp{rQ,t) vs. 5*2(^0, t) for p = 4 
(red open circles) and p = 8 (blue filled circles) with vq/Lz = 0.012, in the middle of the inertial 
range for simulation B. Red, continuous lines represent the intermittent prediction f3p = p — 2C,p 
with C,p given by spatial structure functions; blue dashed lines are the non-intermittent prediction 
/3p = p/3. 

As a further numerical support of ([319]) we now consider temporal behavior of structure 
functions. From ([8]), taking into account the temporal evolution of large scale quantities, we 
expect the temporal scaling Sp{t) ~ t^p with f3p = p — 2(p. With Kolmogorov scaling one 
simply has Pp = p/3 but intermittent corrections are expected to be important, for example 
(3q ~ 2.4 instead of p/3 = 2. Figure [13] shows the scaling of Sp{r,t) vs. S2{r,t) {i.e. in 
the ESS framework) for a particular value of r = Tq = O.OOI2L2. The relative temporal 
exponents Pp/ obtained from the spatial exponents of Fig. [12] fit well the data, while 
non-intermittent relative scaling exponents /3p//32 = p/2 are ruled out. 

The effects of intermittency are particularly important at very small scales. One impor- 
tant example is the statistics of acceleration which has recently been the object of experi- 



mental and numerical investigations 38|, [39|]. For completeness, we briefly recall the main 
results obtained in those studies. 



18 



The acceleration a of a Lagrangian particle transported by the turbulent flow is by def- 
inition given by the r.h.s of ([1]). In the present case of Boussinesq approximation, the 
acceleration has three contributions: pressure gradient, viscous dissipation and buoyancy 
terms. Neglecting intermittency for the moment, dimensional scaling ([MnD implies that 
— Vp ^ uAu ~ i'^^^^((3g6Q)^^'^t^^^ while (3gT ~ (3g6Q. Therefore the buoyancy term in ([T]) 
becomes subleading not only going to small scales but also at later times. Among the other 
two terms, we flnd that, as in standard NS turbulence, the pressure gradient term is by far 
the dominant one, as shown in the inset of Fig. [TH After an initial transient, we have that 
for t > 2r both terms grow with a constant ratio {dzp)rms/{^^ui)rms — 8- 




a/a 



rms 



FIG. 14: Probability density function of the vertical component of the acceleration at time t = It 
(red, bottom tails), t = 2t (green, intermediate tails) and t = 3.8 r (blue, upper tails) normalized 
with rms values. Inset: evolution of arms with time for the three contributions of ([T]): pressure 
gradient dzP (red circles), dissipation I'Aw (green open triangles) and buoyancy term f3gT (blue 
filled triangles). Data from simulation B. 

The inset of Fig. [H] suggests that the temporal growth of arms is faster than t^^^. Again, 
this can be understood as an effect of intermittency which is particularly important at small 







0.86 



39|. 



scales. Indeed, using the multifractal model of intermittency js] one obtains the prediction 



The effect of intermittency on acceleration statistics is evident by looking at the probabil- 
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ity density function. Figure [H] shows that the distribution develops larger tails as turbulence 
intensity, and Reynolds number, increases. This effect is indeed expected, as the shape of 
the acceleration pdf depends on the Reynolds number and therefore no universal form is 
reached. Nevertheless, given the value of R\ as a parameter, the pdf can be predicted again 



using the multifractal model 39 1. 



VI. CONCLUSION 

We have studied spatial and temporal statistics of Rayleigh-Taylor turbulence in three 
dimensions at small Atwood number and at Prandtl number one on the basis of a set of 
high resolution numerical simulations. RT turbulence is a paradigmatic example of non- 
stationary turbulence with a time dependent injection scale. The phenomenological theory 
proposed by Chertkov |3] is based on the notion of adiabaticity where small scales are slaved 
to large ones: the latter are forced by conversion of potential energy into kinetic energy; 
the former undergo a turbulence cascade flowing to smaller scales until molecular viscosity 
becomes important. In this picture, temperature actively forces hydrodynamic degrees of 
freedom at large scales while it behaves like a passive scalar field at small scales where a 
constant kinetic energy flux develops. 

The above scenario suggests comparison of RT turbulence with classical homogeneous, 
isotropic, stationary Navier-Stokes turbulence, in the general framework of the existence of 
universality classes in turbulence. 

By means of accurate direct numerical simulations, we provide numerical evidence in 
favor of the mean-field theory. Moreover, we extend the analysis to higher order statistics 
thus addressing the issue related to intermittency corrections. By measuring scaling expo- 
nents of both velocity and temperature structure functions, we find that indeed they are 
compatible with those obtained in standard turbulence. This result gives further support 
for the universality scenario. 

We also investigate temporal evolution of global quantities, both geometrical (the width 
of mixing layer) and dynamical (the heat flux). The relevant dimensionless quantity in 
RT turbulence are the Rayleigh, Reynolds and Nusselt numbers for which there exists an 



old prediction due to Kraichnan 



19| , known as the "ultimate state of thermal convection" , 



which links the dimensionless number in terms of simple scaling laws. Our set of numerical 
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simulations give again strong evidence for the validity of such scaling in RT turbulence at 
small Atwood number and at Prandtl number one thus confirming how important in thermal 
convection is the role of boundaries which prevent the emergence of the ultimate state. 
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